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ABSTRACT 

A longstanding puzzle of fundamental importance in modern cosmology has been the origin of the 
nearly universal density profiles of dark matter halos found in N-body simulations - the so-called NFW 
profile. We show how this behavior may be understood, simply, by applying adiabatic contraction to 
peaks of Gaussian random fields. We argue that dynamical friction acts to reduce enormously the 
effect of random scatter in the properties of initial peaks, providing a key simplification. We compare 
our model predictions with results of the ultra-high resolution Via Lactea-II N-body simulation, and 
find superb agreement. We show how our model may be used to predict the distribution of halo 
properties like concentration. Our results suggest that many of the basic properties of halo structure 
may be understood using extremely simple physics. 



1. INTRODUCTION 

Hierarchical structure formation is a messy process. 
The assembly of virialized objects in Cold Dark Matter 
(CDM) cosmologies involves countless merger and accre- 
tion events that can convulse the interiors of growing 
halos. N-body simulations have shown that halos can 
have tumultuous lives, in which they are bombarded on 
all sides by infalling clumps of material. 

Out of this seeming chaos, however, emerges remark- 
able regularity. Simulated dark matter halos have fairly 
generic radial profiles, with density scaling as p oc r~ 3 
at large radii, becoming more shallow at smaller radii, 
approaching p ~ r~ x nea r the resolu tion limit of the 
simulations. INavarro et~aT1 (| 19961 119971 hereafter NFW) 
suggested that this behavior is universal among CDM ha- 
los. Subsequent numerical work found qualitatively sim- 
ilar results, producing halos whose density rises steeply 
down to the resolution limits of the simulations (e.g . 
Moore et ai1ll998l: IDiemand et alJl2007l: iGao et al.ll2008t 



Stade l et all 120091: iDiemand et all 120081: INavarro et alj 
20101) . The vast majority of halos show this behavior; 



the exceptions largely appear to correspond to 'bridged' 
halos artificially linked together, or halos that have un- 
dergone recent major mer gers and not yet had time to 
virialize (|Lukic et alll2009l ). 

The origin of this near-universality of halo structure 
has been a longstanding puzzle that has attracted con- 
siderable theoretical attention. Many different types 
of arguments and mechanisms have been advanced to 
account for hal o s' NF W-like behavior. For example, 
iNusser h Shethl (fl99 91 suggested that the shape of the 
halo profile may be related to the shape of the mat- 
ter power spectrum. Other groups have argued that 
tidal disruption of substructure in halos could dynam- 
ically drive the inner profile towards r™ 1 (Syer fc Whitel 



119981: iDekel et all l2003bl lal). While such mechanisms 
could affect halo profile shapes, however, these mecha- 
nisms are not apparently required to produce NFW-like 
profiles. In particular, calculations of monolithic col- 
lapse of halos has shown that the resulting halo s gener- 
ically ha ve NFW-like profiles dHuss et alJll999H More 
recently, iWang fc Whitel (|2009l ) used simulations of Hot 
Dark Matter cosmologies to show that halos forming at 
the cutoff scale of the power spectrum have radial pro- 
files that are fit by the NFW form just as well as CDM 
halos are at comparable resolution. These halos have no 
subhalos in them, and the power spectrum shape on the 
relevant scales is completely different than CDM power 
spectra, and yet the same generic halo profile results. 

Because such a broad class of initial perturbations pro- 
duce collapsed halos with NFW-like profiles, it is useful 
to study particularly simple classes of collapsing halos to 
help identify the importan t mechanisms. Toward s this 
end, in a companion paper (|Lithwick fc Da lal 2010, here- 
after Paper I) we studied the collapse of initially scale- 
free density profiles. This problem admits a similarity 
solution, which allows us to achieve considerably higher 
resolution than is possible using conventional N-body 
techniques. We showed in Paper I that NFW-like profiles 
are not a generic outcome of cold, dissipationless gravita- 
tional collapse in three dimensions. Instead, self-similar 
collapse produces halos with central cusps whose shapes 
depend on the slopes of the initial peaks of the linear 
density. We identified adiabatic contraction as a crucial 
mechanism in setting the overall shape of the halo, and 
showed that the conserved adiabatic invariants governing 
the contraction of mass shells are set during the initial, 
quasi-linear regime preceding collapse into the halo. Us- 
ing this result, we were able to write down a simple toy 
model to predict the final collapsed profiles of halos as 
a function of the initial peaks. In this paper, we will 
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attempt to generalize and apply these results to the for- 
mation of halos in the context of hierarchical structure 
formation. 

This paper is structured as follows. In ^21 we briefly re- 
view the results of our study of self-similar collapse from 
Paper I, and discuss how these results may be generalized 
to peaks that are not scale-invariant. In we compare 
our model predictions with the high-resolution N-body 
simulation Via Lactea-II. We show that our model agrees 
well with the simulation's results, however we find evi- 
dence that dynamical friction significantly modifies the 
structure of the collapsed halo. In 2] we present a sim- 
ple model for how to account for the effect of dynamical 
friction on the statistics of the hierarchy of peaks within 
peaks of Gaussian random fields, again finding excellent 
agreement with the Via Lactea-II simulation. We show 
how our model may be used to predict the distribution 
of halo properties like concentration in SjSJ and conclude 
in gl 

2. HALO COLLAPSE 

Many aspects of halo formation may be understood 
from fairly general considerations of gravitational col- 
lapse, as we will show. We start by reviewing spher- 
ical collapse of peaks. It will prove useful through- 
out the discussion in this section to imagine decompos- 
ing the initial volume around the peak into Lagrangian 
shells, indexed by their Lagrangian radius tl, and to 
build up the collapsed halo profile by summing over La- 
grangian shells. Our discussion and notation here inten- 
tionally parallels that in Paper I, however in this paper 
we will use conventional mass and distance coordinates, 
rather than the self-similar units of Paper I. Besides Pa- 
per I, there have been many previous papers employ- 
ing a similar ap proach toward s understanding halo for - 
mation; see e.g. IPeeblesI (119801); iRvden fc Gunnl (|1987t ); 
IWhite fc Zaritskvl (|1992l ); IDel Popolol (|2009D for exam- 
ples. 

Consider a spherically symmetric peak of the initial 
linear density field <5h n (7L)- Here, S refers to the average 
interior mass overdensity SM/M, not the local overden- 
sity Sp/ p, an d h, is Lagrangian ra dius. The spherical col- 
lapse model (jGunn fc Gottlll972h shows that turnaround 
occurs at the scale where the linear mass overdensity be- 
comes of order unity, that is <5ii n (^c) ~ 1- The linear 
density grows like the linear growth factor, D(a) ~ a at 
early times, so we can easily determine the time at which 
different scales turn around and collapse, given the lin- 
ear density profile. For example, if the linear density has 
some local slope 7, such that locally 5u n oc rr -7 , then 
the expansion factor a when a given scale tl collapses 
behaves as a a tl 1 . Therefore the turnaround radius in 
proper (not comoving) coordinates behaves as r oc r-^ 1+1 . 
From the dependence of the turnaround radius with ini- 
tial Lagrangian radius, we can infer the halo profile under 
various assumptions. 

2.1. Apoapses 

First, for simplicity, let us assume (briefly) that subse- 
quent to turnaround, matter shells remain fixed at some 
constant fraction of the turnaround radius. We refer to 
this as the "frozen" model, since instead of allowing mass 
elements to orbit, we freeze them at their turnaround 
radii. This is essentially identical to the "circular orbit" 
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Fig. 1. — Simplest 1-D models for the halo profile. The dotted 
black line shows the input linear overdensity profile, which was 
adapted from the stacked La grangian profiles o f low-mass halos 
from the N-body simulations of Dalai ct al. ( 2008). It has vanishing 
slope at small radius, and a divergent slope at finite radius due 
to the overdensity becoming negative. The solid curves show the 
predicted profile for various 1-D models. The red curve shows the 
prediction of the frozen model, as time t — ¥ 00. The flat region in 
the linear profile leads to the flat region in the halo profile, while 
the steep region in the linear profile produces the roughly r~ 3 slope 
in the outer halo. The green curve shows the profile if each shell 
lays down interior profiles with minimal tails, without any resulting 
adiabatic contraction. The density is larger than the frozen model 
prediction, but still rolls over in slope quickly. The blue curve 
shows the prediction of the minimal contraction model. Once any 
contraction is taken into account, the roll-over in slope from steep 
to shallow becomes considerably more gradual with radius than the 
frozen model. For comparison, the black dashed curve shows the 
NFW profile, normalized to match the minimal contraction model 
at r_2, the radius where the logarithmic slope dlog p/dlog r = — 2. 
At small radius, the slope of the minimally contracted profile rolls 
over more quickly than the NFW slope. Note that for all the 1-D 
models, the actual outer slope is slightly steeper than -3, since the 
total mass of the halo is finite. 

model of lRvden fc Gunnl (| 19871 ) , who froze mass shells at 
one-half the turnaround radius. The point of this exer- 
cise is to isolate the effect that the initial peak profile has 
on the distribution of orbital apoapses, before shell cross- 
ing occurs. This is important in setting the outer profile, 
though at smaller radii the orbital motions of particles 
lead to significant changes in the mass distribution. 

Making this 'frozen' assumption, we can easily deter- 
mine the density profile using p oc cPr^/cPr. For ex- 
ample, if the linear density is a power-law with slope 7, 
then as noted above r oc tl 1+7 , and the density becomes 
P oc (r-L / r) 2 dr-L / dr oc r~ 9 , where g = 37/(1 + 7) is the 
Fillmore- Goldreich slope described in Paper I. Figure Q] 
illustrates this behavior. In regimes where the slope of 
the linear density is very shallow, 7^0, then the frozen 
model predicts that the final collapsed profile is also shal- 
low, g ~ 0. In contrast, when the initial slope is quite 
steep, 7 — > 00, then the final slope also becomes steep, 
.9^3 (actually slightly steeper when S is not an exact 
power law, i.e. when d"f/dr is nonzero). 

This already explains some of the behavior found in 
halos in cosmological N-body simulations, which tend to 
have steep ~ r~ 3 outer profiles. As can be seen from 
the simple argument above, this is a natural outcome 
whenever the initial linear density peaks acquire steep 
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outer slopes (e.g. when the linear overdensity falls be- 
low zero, a common occurrence among peaks of random 
fields). Put in more physical terms, steep r -3 profiles 
naturally occur whenever the growth of halos slows or 
stops. Indeed, once the halo stops growing, then as 
time progresses the r -3 region simply extends. Essen- 
tially the same argument has previously been used to 
explain the correlations between halo concentration and 
formation time, which underlies the mass-concentration 
relation at low mass jWechsler et all [200l iZhao et al.l 
l2003allbl : ILu et al.ll2006f ). We discuss this point in more 
detail in $5] 

2.2. Periapses and shell profiles 

The frozen model described above assumes that each 
Lagrangian shell lays down a thin annulus of material in 
the collapsed halo, with no shell crossing. This aspect is 
unrealistic, of course - following collapse, material enter- 
ing the halo will execute orbits that cover a finite range 
of radii, and therefore the density profile laid down by 
that material is not a thin annulus, p ^ S(r — r apo ), but 
instead has tails extending to small radii. This affects 
the overall halo profile in multiple ways. First, the slope 
of the interior profile will obviously be steeper than the 
Fillmore- Goldreich slope if the tails from individual shells 
are steep enough. Secondly, because of shell crossing, the 
mass interior to a Lagrangian shell grows following col- 
lapse, which causes contraction and hence steepening of 
the overall profile. 

A full description of these effects requires numerical 
calculation (see Paper I), however we can understand 
some of the basic behavior using simple qualitative ar- 
guments. We can generally expect the density profile 
Psheii(?~) laid down by some shell with total mass Af s h c ii 
to take the form 

Afshoii t(r) 
AsheiiW oc — x - — x P(r pori < r). (1) 

T ^orb 

Each term in this expression is easy to understand. First, 
we expect the density at each scale to behave like mass 
divided by volume, i.e. M^di/r 3 . The middle term cor- 
responds to the fraction of the orbital time that a par- 
ticle spends inside radius r, assuming it reaches radius 
r. Lastly, the third piece is simply the probability that 
the particle is on an orbit whose periapse is inside radius 
r; particles obviously deposit no mass at radii r < r pcr i 
inside of their periapses. 

Next, the time spent at radius r scales as r/v, where v 
is the velocity. Now, as long as the circular velocity pro- 
file is either fiat or declining towards small radius (i.e. 
the total density is p ~ r~ 2 or shallower), then the ve- 
locity of a particle at radii smaller than apoapse will be 
approximately constant, so that t(r) oc r. Therefore, to 
a reasonable approximation, the density laid down by a 
particle behaves as p cx r~ 2 for r pcr i < r < r apo . There 
are small corrections to this near periapse and apoapse, 
but these corrections are unimportant when we sum over 
all particles. Accordingly, the density from a shell be- 
haves like p s heu( r ) ^ r ~ 2 P( r pcii < r), i.e. the density 
profile depends upon the distribution of periapses. 

For example, in spherical collapse where all motion 
is perfectly radial, all particles pass through the origin, 
meaning that P(r pcr i < r) = 1 for all finite r. Then we 



expect the interior profile from each shell to behave as 
r~ 2 , and as long as the Fillmore- Goldreich slope is shal- 
lower than this (i.e. 7 < 2), then the interior profile will 
be r~ 2 . This is indeed the b ehavior found in the full so- 
lutions to spherica l collapse ([Fillmore fc Goldreichlll98'5 
lBertschingerl ll985') . as described in Paper I. 

More generally, however, the orbits of particles in the 
collapsed halo are not perfectly radial, meaning that 
the interior profile is asymptotically shallower than r~ 2 . 
We cannot quantitatively determine the periapse distri- 
bution without characterizing the full orbital structure, 
which requires numerical calculation. However we can 
make useful qualitative arguments. Assuming that the 
density profile is shallower than r~ 2 , then particles at 
small radius r -C r apo have velocities much greater than 
the local circular velocity. We can therefore treat their 
motion as unaccelerated, i.e. straight lines. The distribu- 
tion of the orientations of these straight lines then deter- 
mines the periapse distribution. For example, for random 
orientations with some spread (e.g. a Gaussian distribu- 
tion), the probability to reach radius r simply scales like 
the area subtended by radius r, that is P(r per ; < r) oc r 2 . 
In this case, we would expect that at asymptotically 
small radius, the density profile laid down by individ- 
ual shells should tend towards a consta nt, p s heii — > cons t 
as r — > 0. Note that this disagrees with ILu et all(|2006t) . 
who claimed that an isotropic velocity dispersion leads 
to a r _1 density profile, due to an unphysical ansatz they 
assumed. 

2.3. Shell crossing and contraction 

When material is accreted onto halos, its mass is de- 
posited over a range of radii. This newly deposited mass 
not only adds to the existing density, but also causes 
contraction of the mass already present in the halos. 
To see this, note that the radial action J r = § v r dr oc 
[r x Af(r)] 1 / 2 is an adiabatic invariant for spherically 
symmetric systems. If M at radius r increases due to 
accretion of new matter, then the orbits of particles at 
that radius will shrink in order to hold fixed the product 
r x M(r). This effect is caused entirely by the density 
tails described in the previous subsection; in the absence 
of such tails (i.e., in the absence of shell crossing) the 
mass interior to any shell would remain constant and 
there would be no further contraction following collapse. 

Because the contraction acts to hold fixed rxM,we 
can easily estimate its effects given the r x M profile 
before shell crossing (e.g. at turnaround) and the density 
tails laid down by each shell. The most conservative 
estimate of this contraction arises if we assume minimal 
tails. As noted above, in triaxial potentials, we expect 
shell profiles to behave as p — > const as r — > 0, so the 
weakest tails we expect would be for perfectly constant 
density p = const for all r < r apo . We denote this the 
minimal contraction model. 

In this minimal contraction model, a Lagrangian mass 
shell Ml, of width dM^ and apoapse r apo (Mi J ) lays down 
a mass profile M(r) = dMz,f{r/r apo ), where 

/(r/r apo ) = j fc) ' r<r a P o , (2) 
I 1, r > r apo 

We compute the total mass profile by summing over all 
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shells. This requires knowing r a p Q (ML), which we can 
estimate by assuming that the product r x M is constant. 
At turnaround, this product is 

F(M h ) = M L x r ta (M L ) cx A^ /3 /4„(M L ), (3) 

where we have used r ta oc r^/S^. 

Let us write r apo (ML) as the apoapse for shell Ml, and 
the inverse function M a (r) which gives the Lagrangian 
shell Ml whose r apo = r. Then we can easily write down 
the total (contract ed) mass profile at t = oo as a sum 
over all shells (e.g. iRvden fc Grmr][l987t ILu et apjOOg) 

Inserting our expression for the shell profile, Eqn. @, 
gives 

r°° / \ ^ 

M(r) = M a (r) + / dM L — - . (5) 

Differentiating this expression with respect to r gives 

dM 3 r , w . 
— = - M r - M a r 
ar r 

= -[M - F-\Mr)}, (6) 
r 

where F^ 1 is the inverse function of the expression in 
Eqn. This is a simple ordinary differential equation 
for M(r) which may easily be solved for any linear den- 
sity profile c>ii„(r L ). 

As the name implies, the minimal contraction model 
provides a lower limit on the effects of adiabatic con- 
traction. If the shell profiles have more mass at small 
radius than this model assumes, then contraction can be 
considerably stronger. To illustrate this, we provide an- 
other example toy model, (very) loosely motivated by 
our earlier study of self-similar triaxial collapse. We 
found that in certain cases, the shell profiles had inner 
slopes d(log/9 sh cii)Mlogr) ~ ±d(log p tot )/d(\ogr). If we 
assume that all shells have such profiles inside of their 
apoapses, then following the same reasoning that led us 
to Eqn. ©J it is straightforward to show that the mass 
profile satisfies 

dM _ 3M M - F^ 1 (M r) 

~~dT ~ r M + F~ 1 (M r) ' ( ' 

We label this as the p 1 / 2 model. We stress that both 
of these toy models are not meant to be considered as 
serious descriptions of the shell profiles; rather they are 
meant to be illustrative, since they allow the halo mass 
distribution to be computed by solving simple ordinary 
differential equations. 

Figure [T] illustrates the effects of contraction. The pre- 
dicted halo profile is significantly modified compared to 
the frozen model, even with minimal contraction. As the 
figure shows, the transition from steep to shallow slopes 
becomes much more gradual, once the effects of contrac- 
tion are taken into account. Even though the input peak 
profile sharply transitions to a flat, 5 ~ const behavior, 
similar to a top-hat perturbation, the adiabatically con- 
tracted halo profile is very similar in form to NFW, even 
with minimal contraction. It is therefore not entirely 



surprising that NFW-like profiles arise in contexts like 
monolithic collapse or simulations with truncated power 
spectra. 

It is worth noting, however, that our models do not 
generically predict power-law central cusps in the final 
halo profile p(r). Only for power-law initial profiles, 
S oc r -7 do we find power-law cusps, p oc r~ 9 (see 
Paper I). More generally, for peaks with 5 — > const 
as tl — > 0, our models predict only logarithmic diver- 
gence of the halo density p as r — > 0. The logarith- 
mic slope dlogp/dlogr — » as r approaches 0, however 
the roll-over in slope occurs extremely slowly over many 
dec ades in radius. Recent high resolution simulations 
(e.g iStadel et all 120091 : iNavarro et all 12010( 1 have found 
similar behavior, in the sense that their halo profiles ap- 
pear be tter described with r olling profiles like the Einasto 
profile (|Merritt et a l. 2005) instead of power-law cusps. 
We return to this topic in $6l 

To summarize this section, we have described a sim- 
ple method to translate from the initial peak to the final 
halo, essentially by adiabatically contracting the linear 
density profile and adopting a prescription for the dis- 
tribution of orbits. In the next section, we compare this 
model to results of a high resolution cosmological N-body 
simulation. 

3. COMPARISON WITH N-BODY SIMULATIONS 

The model discussed in the previous section may seem 
overly simplistic. It employs arguments based on spheri- 
cal symmetry, and describes halo assembly as an orderly 
process, resulting in an effectively stratified structure in 
which the orbits of particles within the halo reflect the 
locations where those particles originated. In contrast, 
the assembly of halos observed in cosmological simula- 
tions is often violent, involving discrete, stochastic ac- 
cretion events that frequently take the form of major 
mergers. For this reason, it is unclear whether a model 
that assumes that all halo material evolves adiabatically 
is capable of describing the messy reality of structure 
formation. 

Previous work does, however, give us reason to hope 
that our simple model may nonetheless be useful. Results 
from multiple simulations have shown that hierarchical 
assembly, including major mergers, does not completely 
obliterate all pre - existing structure within halos (e.g. 
iKazantzidis et all 120061: iValluri et all I2007D . More re- 
cently, iWang et al.l (|2010( 1 analyzed high-resolution sim- 
ulations of individual halos, and showed that their halo 
structure is indeed somewhat stratified, with a clear gra- 
dient of accretion time with radius. On the other hand, 
they also showed that major mergers can disrupt this 
stratification, by bringing in fresh material into the halo 
core. Therefore, it is not clear whether our model, with 
its simplistic assumptions, would be relevant for realistic 
halos. 

To test the assumptions and predictions of our simple 
model, we have analyzed results of the ultra- high res- 
olutio n Via Lactea-II (VL2) simulation ([Diemand et al.l 
2008). This simulation resolves the Lagrangian region 
of a typical Milky Way dark matter halo with just over 
one billion particles of mass 4,100 M Q , and follows its 
evolution and formation in a cosmological environment 
(40 Mpc), most of which is covered only with lower res- 
olution (higher mass) particles. The simulation begins 
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Fig. 2. — Initial linear density profile of the VL2 peak. The thick 
solid black curve shows the interior overdensity profile of the La- 
grangian volume at z = 104.3, linearly evolved to redshift z = 0. 
Lagrangian radius is measured relative to the Lagrangian centroid 
of particles in the halo's main progenitor at z = 17.88. The thick 
dashed black curve shows the average profile about the Lagrangian 
centroid of the particles inside r200 at 2 = 0. For comparison, 
the upper (blue) line and shaded region indicate the mean and 
dispersion in the highest subpeak profile expected from supremum 
statistics (see ij4jl, while the lower (gray) line and shaded region 
depicts the mean and dispersion expected for the peak profile ne- 
glecting subpeaks (BBKS). 

at redshift z = 104.3, and outputs 400 snapshots in 
time ending at z = 0. The snapshots were processed 
to determine halo properties including mass, centroid, 
and mean velocity, allowing us to reconstruct the halo's 
detailed assembly history. At z = the main halo 
has r2oo = 402 kpc (the radius enclosing a density of 
200/9M = 200f2 M /O cri t) and M 200 = 1-92 x 10 12 M Q . Its 
density profile is well fit by an Einasto profile with pa- 
rameters p_ 2 = 9.91 x 10 S M Q kpc~ 3 , r_ 2 = 28.9 kpc, 
and a = 0.142. Assuming that this simulation provides 
a typical example of halo formation for galaxies like the 
Milky Way, its high resolution should allow us to test the 
key assumptions and predictions of our model. 

The most important assumption of our model that we 
would like to test is adiabaticity of the orbits. A di- 
rect check of this for over 1 billion particles with 400 
snapshots would be a computational challenge, which we 
defer to future work. Instead, we will attempt an easier 
exercise which worked well in our study of self-similar 
triaxial collapse in Paper I. As mentioned in Sj^l we ex- 
pect that the product r x M(r) should be approximately 
conserved if the radial action J r = <f v r dr is an adia- 
batic invariant. In Paper I, we computed the average of 
this quantity, f x M{r) for entire shells of particles, and 
showed that the shell average is not only well conserved, 
but furthermore may be predicted from the linear density 
profile of the initial peak. Because this is much simpler 
to compute, we perform a similar test for the VL2 halo. 

The first ingredient needed is the linear density profile 
of the initial peak that collapses to form this halo. This 
is plotted as the solid black curve in Figure [2j Given this 
peak profile, we predict each shell's invariant r x M us- 
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Fig. 3. — Adiabatic contraction. The solid black curve shows 
the median radius r at z = for particles in narrow Lagrangian 
shells, as a function of Lagrangian radius rx. The dashed blue 
curve shows the predicted final radius for each shell, under the 
assumption that each shell adiabatically contracts, with adiabatic 
invariants determined from the linear density profile of the initial 
peak at z = 104.3. For comparison, the dotted red curve shows the 
predicted final radius in the frozen model, in which r oc r ta with 
no subsequent adiabatic contraction. 

ing the spherical collapse model. At early times, prior 
to shell crossing, the mass enclosed within each shell 
is a constant, = (47r/3)prL 3 , and the spherical col- 
lapse model predicts the shell's turnaround radius to be 
r ta = 0.6rL/5(rL). We assume that each shell's adiabatic 
invariant is equal to r ta x Mi,. 

Given each shell's adiabatic invariant, we can compute 
the shell's average final radius within the halo at z = 0, 
given the halo's mass profile M(r): we simply find the 
radius where rxM(r) is equal to r ta x Ml for each shell. 
This is plotted as the dashed blue curve in Figure [3l For 
comparison, the solid black curve in Fig. |3] shows the 
median radius r measured at z = for each shell. The 
level of agreement between the two curves is striking, 
especially since neither has any freedom to be adjusted. 

It is important to stress that this test is not circular: 
this is a direct check that shells conserve r x M(r), the 
prediction of adiabatic evolution. For a given M(r) pro- 
file, there are many ways that the Lagrangian shells could 
add up to give the total mass. For example, if the halo 
had violently relaxed during formation, such that all the 
orbital actions became randomized, then the final orbital 
actions (and hence orbital radii) would be unrelated to 
original Lagrangian location, and we would expect the 
median r to be independent of r^. This clearly does not 
occur within the VL2 halo. 

To illustrate that this test could have failed, we plot 
in Figure 2] how this comparison would have looked if 
we had used a different Lagrangian centroid. The <J(7"l) 
profile plotted in Fig. [2] was measured relative to the La- 
grangian centroid of the progenitor halo at high redshift, 
z = 17.88. If we instead use the Lagrangian centroid of 
all the particles in the halo at z = 0, we obtain the result 
shown in Fig. @] Here, the agreement between the pre- 
dicted and measured radii at z = is reasonable at large 
radii, but much worse near r = 0. There are two reasons 



6 



0.01 



0.001 




Fig. 4. — Like figurcO but using the 2 = Lagrangian centroid. 
The disagreement at small r shows that the mass at small radius, 
r ~ 0, does not originate near the centroid of the overall Lagrangian 
volume, but instead originates near the subpeak that was used in 

Fig. m 

for this. First, the linear density profile S(ri,) changes 
when we measure Lagrangian radii relative to a different 
centroid, as illustrated by the difference between the solid 
and dashed curves in Fig.[5J More importantly, however, 
the particles originating near the z = Lagrangian cen- 
troid are not the particles that are found near r = in 
the final collapsed halo. Rather, the mass at small r in 
the final halo largely originated near a sub-peak within 
the total Lagrangian volume, that gave rise to the main 
progenitor halo at high redshift. This sub-peak is off- 
center within the overall Lagrangian volume (see Fig. [5]), 
but the mass within this off-center peak eventually falls 
towards r = th rough processes like dynamical friction . 
Earlier work (e.g. iDiemand et al.l2005HWang et al.l2010h 
has already shown that the mass at small radius within 
halos typically collapses at high redshift. Our result is 
potentially much more powerful: we now have a means 
of quantitatively calculating where material will occur 
within the halo, as a function of time. 

The level of agreement shown in Fig.[3]is even more re- 
markable when we consider how simplistic the prediction 
is. To reiterate^ we take the spherically averaged linear 
density profile S(r^), apply the spherical collapse model 
to predict each shell's turnaround radius rt a) and then as- 
sume that the product r ta x Ml is an adiabatic invariant. 
The formation of the VL2 halo is highly nonsphericaQ, 
as is typical in CDM cosmologies, so one potential area 
of future work could be to explore whether the adiabatic 
invariants may be predicted more precisely by accounting 
for triaxiality. 

This simplistic, spherical approach is already adequate 
to make considerable progress towards understanding 
and predicting the mass distribution within the collapsed 
halo. If we know how to calculate the typical location of 
each shell within the halo, then we need only add a pre- 
scription for the individual shell profiles to predict the 

1 The nonspherical collapse of this halo may be seen in the movie 
http : / /astro .berkeley . edu/~mqk/VL2/movie_10M_withf ly . mp4 




Fig. 5. — Halo mass profile. The solid black curve is M(r) mea- 
sured from VL2, the dashed blue curve is our p 1 / 2 toy model Eqn. 
ff). and the dotted red curve is the minimal contraction model 
Eqn. J6j- 

total mass profile. In $2j we described two examples 
of toy models for the profiles of shells: a minimal con- 
traction model, in which each shell deposits density as a 
step function, /9 s hcii oc 0(r apo — r), and a slightly more 
elaborate model in which each shell deposits mass with 

1 /2 

a profile p s i lc u oc p tot Q(r apo — r). These are both obvi- 
ously very crude treatments and are no replacement for 
a detailed understanding of the orbital structure within 
halos; we employ them here simply because they allow 
us to solve for the total mass distribution using ordinary 
differential equations. We plot the predictions of these 
toy models in Figures [5] and |6l Note that, unlike the plots 
in Figs. |3] and HI this calculation does not use the mea- 
sured M(r) profile from the VL2 halo at z = 0. Rather, 
these toy models solve self-consistently for the total mass 
distribution, placing each shell at the required radius r 
given each shell's adiabatic invariant and the mass pro- 
file obtained by summing over all shells. Fig. |6] shows 
that this toy model is only a very crude approximation 
to the individual shell profiles, but Fig. [5] illustrates that 
this is already adequate to predict the halo's radial mass 
distribution reasonably well. 

Before closing this section, we stress that VL2 is just 
one halo. The remarkable agreement between our model 
predictions and this simulation's results could, in princi- 
ple, be a fluke, if (for example) this halo had an atypi- 
cal forma tion history. Further study of additiona l exam- 
ples (e.g. iStadel et alJ 12009b iNavarro et all 12010( 1 would 
be useful to verify the generality of our results. 

4. PEAKS AND SUB-PEAKS OF GAUSSIAN RANDOM 
FIELDS 

In the previous section, we showed that the radial 
structure of the VL2 halo may be understood quite sim- 
ply. Lagrangian shells evolve over time in a manner that, 
on average, conserves their radial actions. These adi- 
abatic invariants may be predicted in a simple way by 
applying the spherical collapse model to the spherically 
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Fig. 6. — Mass profiles for two particular Lagrangian shells. 
The solid curves show the profiles measured from VL2, the dot- 
ted curves are the minimal contraction toy model, and the dashed 
curves are the p 1 / 2 toy model. Neither toy model captures the shell 
profiles well, especially at small radius, however the p 1 / 2 model is 
a close enough approximation to reproduce the overall mass profile 
(see Fig. [5j. 

averaged profile of the peak that collapses to form this 
halo. Given this result, even a crude treatment of the 
mass profiles deposited by Lagrangian shells suffices to 
describe the overall mass profile within the halo reason- 
ably well. Consistent with previous work, we also found 
evidence for processes like dynamical friction that can 
transport material within early-collapsing subregions to- 
wards the halo center. 

Given this success in understanding one particular 
halo, what can we say about halo structure more gener- 
ally? Our key result is that the profiles of collapsed halos 
may easily be understood in terms of the properties of 
their precursor peaks. For most cosmologies of interest, 
the properties of th e initial peaks are w ell described by 
Gaussian statistics (B ardeen et al.lll986l ). 

In our model, the radial profile p(r) within the halo 
is determined by the linear overdensity profile S(r^) 
of the initial peak. It is straightforward to determine 
the statistics of the interior <5(riJ profile for Gaussian 
random fields. For example, suppose that we have a 
peak on scale r p k that collapses to make a halo of 
mass M ~ (47r/3)prp k , and suppose that the peak has 
height <5(r p k) = 5 p k and derivative dE / dr^r-p^) — <5 pk 
on this scale. We would like to compute the condi- 
tional probability distribution for the density at interior 
radji given the peak height and slope at the outer scale, 
P(S( r L)\S P k, S' pk ). In general, the conditional probabil- 
ity distribution P(X\Y) for Gaussian variables X and 
Y (with zero mean) is also a Gaussian, with mean and 
variance 

(X\Y) = (XY)(YY)- 1 Y (8) 
a 2 xlY = (XX)-(XY)(YY)- l (YX). (9) 

In our case, X would correspond to the small-scale den- 



sity <5(rL), and Y would correspond to the density <5 p k 
and slope S' pk on scale r p k- Given a power spectrum, the 
required covariances are easy to compute. In Fig. [2l the 
dotted line and gray dashed band illustrate the mean and 
dispersion of the 5(rt,) profile when we condition on the 
value and slope of 8 at the top-hat radius. The mean pro- 
file quickly plateaus inside the top- hat radius < tth, 
where M200 = (47r/3)pr^ H . In addition, there is consid- 
erable scatter about the mean profile, increasing towards 
smaller radii. _If we use the z = Lagrangian centroid, 
the measured S profile (dashed line in the Figure) is quite 
consistent with this BBKS prediction. Previous workers 
have made similar arguments for t he expected shape of 
the Ga ussian peaks that form halos. Floffman & Shaham 
(|1985| ) assumed that peak profiles follow the unsmoothed 
matter correlation function, whic h BBKS later sh o wed is 
not typical for G aussian peaks. iRvden fc Gunnl ()1987f ) 
and IRvden! (|1988[ ) used the average BBKS profile, how- 
ever they used a smoothing scale much smaller than the 
top-hat scale of the halos, which gives a profile extremely 
atypical for the pea k s tha t collapse into halos. More 
recently, iDel Popolol (|2009| ) employed the mean BBKS 
profile to describe the typical profiles of Gaussian peaks. 

As we have seen in $31 however, the linear profile cen- 
tered on the overall Lagrangian centroid may not be rel- 
evant for computing the structure of the collapsed halo, 
because processes like dynamical friction can drag off- 
center subpeaks towards the halo center. We found good 
agreement between the measured and predicted adiabatic 
invariants if we instead used the Lagrangian position of 
the sub-peak that formed the earliest halo progenitor at 
z = 17.88. As Fig. [2] shows, the profile centered on this 
subpeak is quite different than the mean BBKS profile 
used in earlier works. In general, we would expect sim- 
ilar behavior in other halos as well. The hierarchy of 
peaks within peaks expected for CDM cosmologies sig- 
nificantly modifies the structure that we would naively 
calculate using the mean peak profile. 

One simple way to account for this effect is simply to 
graft the profile of the highest subpeak onto the overall 
peak profile. So if we wish to compute the density <5(rL) 
for some Tl that is much smaller than the peak size, in- 
stead of using the value of S(ri,) centered on the halo's 
Lagrangian centroid, we instead find the largest S for all 
the sub-volumes of size r su b = within the halo's La- 
grangian volume. This is effectively what we have done 
in Fig. [2] by centering on the z = 17.88 sub-peak. Of 
course, this sub-peak will have subsubpeaks inside of it, 
but with a hierarchy of such grafts, we can construct the 
effective peak profile that accounts for effects like dynam- 
ical friction. 

The statistics of the highest subpeaks are elementary 
to compute, as we illustrate with a simple example. Sup- 
pose that a; is a Gaussian random variable with zero mean 
and unit variance, of which we have N independent sam- 
ples. The probability that y exceeds any one of these x 
samples is 

ft W-jCs* _1 -?*(^)' <10) 

and so the probability that y exceeds all N of the x values 




N 

Fig. 7. — Mean and dispersion of the supremum, x sup , the largest 
value of N random Gaussian numbers with zero mean and unit 
variance. 
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Pn(v) = [Pi(y)] N « exp (-Ny- l e-^' 2 /^\ , (11 



where the last approximation is valid in the limit of large 
N . The double-exponential form of the extreme-value 
probability is not specific only to our Gaussian distribu- 
tion, but in fact arises generica lly for any parent distribu - 
tion that is sufficiently steep (jBhavsar fc Barrowlll985l ). 
The differential probability that y is the supremum of the 
N samples is dP^/dy. In Figure [3 we plot the typical 
supremum as a function of sample size N . As might be 
expected for a Gaussian parent distribution, the largest 
value grows roughly like (log N) 1 / 2 in the limit of large 
N. 

It is straightforward to apply this expression to our 
problem of computing the distribution of the highest 
subpeaks. Here, Pi is simply the BBKS probability de- 
scribed above, and N = (r p k/r su b) 3 is the number of 
independent samples of size r su b inside a peak of size 
r pk- By construction, we account for the correlations 
between the density on scale r su b and the density and 
derivatives on scale r p k, but we neglect any additional 
spatial correlations among the subpeaks beyond this. 

We caution that, in general, this expression will over- 
estimate the height of the effective S(r^) profile, and will 
underestimate its scatter. This is because dynamical fric- 
tion is not always effective, especially at low subhalo 
mass - which is why so much substr ucture persists in 
CDM halos (e.g. iDiemand et"all 12008ft . We are implic- 
itly assuming that the highest subpeak of size r su b will 
always find its way to the origin, but for M su b -C A/haio, 
the dynamical friction timescale may exceed the halo's 
lifetime, meaning that this sub-peak would have to be 
carried to r = by some larger scale structure. Clearly, 
it is not always the case that, for example, the high- 
est peak on scale r = 10 -3 arises inside of the highest 
peak on scale r = 10~ 2 . A more careful treatment would 



be considerably more complicated, however, so we have 
opted to make this simplifying assumption. 

Bearing this caveat in mind, the blue curve and shaded 
area m Fig. [2] plots the mean value and dispersion of the 
highest sub-peak, calculated using supremum statistics. 
As with the BBKS curve, we have conditioned the prob- 
ability on the value and derivative of the linear overden- 
sity 5 on the top-hat scale corresponding to the measured 
mass M2oo- The agreement between the predicted curve 
and the actual measured peak profile over decades in ra- 
dius is striking. 

Another remarkable feature is that the scatter in the 
effective profile is much smaller than the scatter in the 
BBKS profile. This occurs because the the dispersion in 
the largest value of a sample is considerably smaller than 
the dispersion of the sample as a whole, in the limit of 
large N. Note that this level of dispersion is not indica- 
tive of the full scatter in peak profiles for halos of a fixed 
mass. In this example, we have constrained the profiles 
to match the value and slope of the VL2 peak profile at 
the largest plotted radius. Other halos with comparable 
mass will correspond to peaks with different heights and 
slopes at the boundary, and hence their internal profiles 
will show different behavior. 

To illustrate this, we plot in Fig. [5] the predicted pro- 
files for initial peaks that have outer slopes different than 
the VL2 peak. The blue and red solid curves show the 
predicted halo profiles for peaks with outer slopes either 
2x larger, or 3x smaller, than the VL2 peak. For com- 
parison, the dotted curve shows the NFW profile, 
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and the dashed curves show the Einasto profile 
(jMerritt et al.ll2005D . 
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for a = 0.13 and 0.1 7, covering the rang e found in the 
Aquarius simulations (|Navarro et al.ll2010f) . As the figure 
shows, varying the outer profile slope of the initial peak 
over a reasonable range spans the profile shapes found in 
high resolution simulations. 

5. CONCENTRATIONS 

The previous sections described a simple model for 
the structure of cosmological halos, combining Gaussian 
statistics with prescriptions for physical effects like dy- 
namical friction and adiabatic contraction. Such a model 
can have many obvious applications. We illustrate with 
one example in this section. 

In our model, the shape of the final collapsed halo pro- 
file is simply related to the shape of the linear density 
profile of the original peak. One of the important param- 
eters used to describe halo profiles is the concentration, 
which we can define as the ratio between the halo's radius 
?™20u j an d r_2, the radius where the local density slope 
is dlogp/dlogr = —2. Many papers have attempted 
to quantify t he typical halo concentratio n as a function 
of mass (e.g. iMuhoz-Cuartas et al.ll2010L and references 
therein). Using our model, we can try to predict this 
c(M) relation. 
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Fig. 8. — Density profile for peaks with the same height as the 
VL2 peak, but with different outer slopes (solid red and blue 
curves), using the p 1 / 2 toy model. For comparison, the dotted 
black curve shows the NFW profile, while the two dashed black 
curves are the Eina sto profile for a = .13 and 0.17, the range of 
the Aquarius halos (Navarro et aL 2 0101) . 

The halo concentration basically measures how quickly 
the halo slope rolls over from near -3 to -2, and in our 
model this is controlled by the outer height and slope of 
the initial peak (mainly the latter). The distribution of 
the heights of the peaks producin g halos is set by the dis- 
tribution of collapse thresholds. iBond fc Mversl (|1996l ) 
proposed a simple ellipsoidal collapse mod el to predict 
the di stribution of collapse thresholds, a nd iSheth et al. 
(2001) found that the predictions of the IBond fc Mvers 
model are reasonably approximated by the fitting func- 
tion 

1+13' " 



(14) 



where 5 C — 1.686 is the spherical collapse threshold, 
a 2 (M) is the variance of linear density fluctuations 
smoothed on mass scale M, and the fit ting function pa- 
ramet ers are j3 = 0.47 and 7 = 0.615. iRobertson et al.1 
(2009) showed that this expression is in broad agree- 
ment with the heights of peaks that produce halos in 
ACDM simulations, so we will adopt it here. More im- 
portantly, we also need the typical outer slopes of the 
same peaks. This has been less well quantified in the 
literature. For high mass halos, with a(M) <C S c , we can 
use simple Gaussian statistics to predict the outer slopes 
5' = dS/dlogrt,. The typical slope for a peak of height S 



is 



(6'\6) 



P(S'\5)dS', 



(15) 



where P is the conditional Gaussian probability distri- 
bution. Note that this expression is different than Eqn. 
©, since we require that the slope 5' < 0, since other- 
wise this peak would not collapse on mass sc ale M but 
instead some larger mass (jBardeen et al.lll986l ). This ex- 
pression works well for high mass halos, M ^> M±, where 
M+ is the characteristic nonlinear mass scale satisfying 
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Fig. 9. — Mean concentration relation c(M). The solid black 
curve shows our model (see text), while the dotted red curve 
shows the prediction from using the naive peak slopes from 
Gaussian statistics. For comparison, the dashed blue curve 
shows the power law c oc M ~ 097 found in N-body simulations 
(Mu noz^uartas et £01120101 ). 

a(M+) = S c . However, this significantly underestimates 
the magnitude of the outer slopes for low mass halos, 
M < M*, which are considerably steeper than Eqn. (|15|) 
predicts. This appears to occur for the same reason that 
these low-mass halos are anti-biased: such halos tend 
to avoid high-density regions. More precisely, the peaks 
that produce low mass halos tend to occur within under- 
dense environments, presumably because similar peaks 
inside of overdense regions do not lead to low mass ha- 
los, but rather are incorporated into higher mass halos. 
This is the reason why the average peak profile plotted 
as the dotted curve in Fig. Q] changes sign. As the fig- 
ure shows, the outer profiles for such halos are extremely 
steep, which we argued in i j2.1l was the origin of the high 
concentrations of these halos. 

Given the lack of a theory to describe this effect, we 
have attempted to measure the typical outer slopes of 
initial peaks. We performed a low-resolution ACDM sim- 
ulation for a volume 256 hr 1 Mpc on a side, identified 
halos at z — down to mass M 1O 12 /i _1 M0, and 
stacked the linear density profiles of the halos' precursor 
peaks in narrow mass bins. Very crudely, we found that 
a rough scaling 
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reasonably described the mass range that we measured. 

Given this expression, and our expression for the av- 
erage peak height as a function of mass, we use Eqn. 
( TTj) to predict the average peak profile, and then Eqn. 
1 7]) to predict the halo profile, from which we measure 
the concentration. Figure [9] shows the result. The solid 
black curve in the figure is our prediction, while the 
dashed blue curve shows the sca ling c oc M~° 097 , as 
found in cosmological simulations (|Mufioz-Cuartas et al.l 
l20ll . For comparison, the red dotted curve shows how 
the concentration would depend on mass if the outer 
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slopes followed simple Gaussian statistics. The agree- 
ment between our prediction and the results of simula- 
tions is, once again, quite good. At low masses, this 
model appears to underpredict the concentrations, but 
this is likely because Eqn. (fl6|) underestimates the slopes 
of halos below M < 10 12 M Q . At high masses, note that 
our concentrations appear to saturate near values c ~ 4. 
This occurs because, at high masses, peaks no longer are 
anti-biased, and hence their outer slopes are given by the 
shallow values predicted from Gaussian statistics. 

6. CONCLUSIONS 

In this paper, we have presented a model to explain 
the origin of the nearly universal density profiles of dark 
matter halos found in N-body simulations. We argued 
that the physics setting the shape of halo radial profiles 
is extremely simple. We find that adiabatic contraction 
sets the basic shape of the halo profile, and that the con- 
served quantities in this contraction, i.e. the adiabatic 
invariants, are determined by the linear density profile 
of the initial peak. We further argue that, because of 
dynamical friction, the hierarchy of peaks within peaks 
significantly alters and simplifies the effective peak pro- 
file setting the adiabatic invariants and hence the halo 
profile. 

We have compared our model predictions to N-body 
simulations, and found striking agreement. In particu- 
lar, the detailed mass distribution of the high-resolution 
Via Lactea-II halo is quite consistent with our model, 
and provides strong evidence for the importance of both 
adiabatic contraction and dynamical friction. We then 
showed how this model may be used to predict the statis- 
tics of halo properties, focusing on the example of the 
mean concentration relation c(M). 

Our model, if correct, could have additional impli- 
cations beyond what we have discussed so far. One 
example that has attracted considerable attention is 
the asymptotic inner profile of dark matter halos, 
whic h can be import ant for dark matter ann i hilation 
(e.g. iBergstrom et al.lll998t iKuhlen et all 120091 : iKuhlenl 
[20iq iReed et al l I2010D . For CDM-like power spectra, 



our models do not produce power-law central cusps, 
but instead diverge logarithmically. The local slope 
d log p/d log r rolls over very slowly with radius, asymp- 
totically approaching zero slope at r = 0. As we 
noted above, recent high resolution simulations have 
found qualitatively similar behavior, where the halo slope 
d\ogp/d\ogr rolls smoothly down to the resolution lim- 
its of the simulations. As we showed in Figure [U the 
Einasto profiles used to fit these halos are quite similar 
to the predictions of our model. 

As we have emphasized in this paper, our model is not 
complete, because it lacks a detailed understanding of the 
mass profiles deposited by Lagrangian shells. For sim- 
plicity, we have adopted a simple ansatz that is loosely 
motivated by our previous study of self-similar triaxial 
collapse (Paper I), however in the future we intend to 
examine the orbital distribution within halos in more de- 
tail. Our study of self-similar collapse has already shown 
that the evolution of orbits in time- varying triaxial po- 
tentials is quite interesting and can often be surprising, 
for example in the different ways that box orbits and loop 
orbits respond to changes in the potential. This is the 
subject of ongoing work. 

Lastly, we note that it would be useful to extend the 
comparison of our models to other high resolution halo 
simulations. In particular, it would be worthwhile to 
explore the limits of this model, which so far appears to 
work surprisingly well in regimes (e.g. r ~ r v j r ) where we 
might naively have expected it to fail. One regime where 
this model is likely to break down is the case of nearly 
1:1 major mergers in halos, in which both progenitors 
have comparable central density and both contribute to 
the central core. It would be interesting to see whether 
we can find adiabatic invariants to describe the dynamics 
even in such extreme cases. 
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